-module(opt).
-export([min1_list/2, minn_list/2, min1/1, min1/3, minn/1, minn/3, minn_man/2, minn_para/2, minn_sx/2, 
         max1_list/2, maxn_list/2, max1/1, max1/3, maxn/1, maxn/3, maxn_man/2, maxn_para/2, maxn_sx/2, 
         test/0, perf_test/0 ]).
-import(make_fun).
-define(ITERATIONS, 10000).
-define(EPS, 0.0001).


%%%%% F/1 in list

% in: F - target function, list of values for which F should be calculated
% out: value from given list, for which F result is the smallest
min1_list(_, [H | []]) ->
    H;
min1_list(F, [H | [TH| T]]) ->
    FH = F(H),
    FTH = F(TH),
    if 
        FH < FTH ->
            min1_list(F, [H | T]);
        true ->
            min1_list(F, [TH | T])
    end.


%%%%% F/n in list

% in: F - target function, list of points for which F should be calculated
% out: points from given list, for which F result is the smallest
minn_list(_, [H | []]) ->
    H;
minn_list(F, [H | [TH| T]]) ->
    FH = apply(F, H),
    FTH = apply(F, TH),
    if 
        FH < FTH ->
            minn_list(F, [H | T]);
        true ->
            minn_list(F, [TH | T])
    end.
    

%%%%% on R

% in: F - target function; iteration count, precision if neccessary
% out: {confirmed, X} if local minimum found in given iterations, {presumably, X} if not
min1( F ) ->
    min1_( F, 0.0, ?ITERATIONS, ?EPS). 

min1( F, ITERATIONS, EPS ) ->
    min1_( F, 0.0, ITERATIONS, EPS). 


min1_( _, X, 0, _) ->
    {presumably, X};
min1_( F, X, I, Eps) ->
    Now = F( X ),
    Left = F( X - Eps ),
    Right = F( X + Eps ),
    if
        Left < Now ->
            min1_( F, X + min1_e_( F, X, -Eps, Left), I-1, Eps);
        Right < Now ->
            min1_( F, X + min1_e_( F, X, Eps, Right), I-1, Eps);
        true ->
            {confirmed, X}
    end.

min1_e_( F, X, E, Last ) ->
    New = F( X + 2*E ),
    case New < Last of
        true ->
            min1_e_( F, X, E*2, New );
        false ->
            E
    end.


%%%%% on Rn

sum_of_d_square(L1, L2) ->
    lists:sum(lists:map( fun({A, B}) -> (A-B)*(A-B) end, lists:zip(L1, L2) )).

x_to_list( X, L, I ) ->
    lists:sublist(L, I-1) ++ [ X ] ++ lists:nthtail(I, L).

% in: F - target function; iteration count, precision if neccessary
% out: {confirmed, [X1, X2, ...]} if local minimum found in given iterations, {presumably, [X1, X2, ...]} if not
minn( F ) ->
    N = proplists:get_value( arity, erlang:fun_info( F )),
    Xs = [0 || _ <- lists:seq(1, N)],
    minn_search_( F, Xs, ?ITERATIONS, ?EPS).
minn( F, ITERATIONS, EPS ) ->
    N = proplists:get_value( arity, erlang:fun_info( F )),
    Xs = [0 || _ <- lists:seq(1, N)],
    minn_search_( F, Xs, ITERATIONS, EPS).

minn_search_( _, Xs, 0, _) ->
    {presumably, Xs};
minn_search_( F, Xs, ITERATIONS, EPS) ->
    N = proplists:get_value( arity, erlang:fun_info( F )),
    Xs_new = minn_inner_loop_(F, N, Xs, [], 1),
    SODS = sum_of_d_square( Xs, Xs_new ),
    if  SODS < EPS -> 
            {confirmed, Xs_new};
        true -> 
            minn_search_( F, Xs_new, ITERATIONS, EPS)
    end.

replace_first(From, To) ->
    From ++ lists:nthtail(length(From), To).
    
minn_inner_loop_(F, N, Xs, Xs_new, I) ->    
    if 
        I == N+1 -> 
            Xs_new;
        true ->
            Fi_F = fun (F_, I_, Xs_) -> 
                        fun ( X ) -> 
                                apply(F_, x_to_list(X, Xs_, I_) )
                        end 
                 end,
            Fi = Fi_F( F, I, replace_first( Xs_new, Xs) ),
            {_, Xi} = min1( Fi ),
            minn_inner_loop_(F, N, Xs, Xs_new ++ [Xi], I+1)
    end.


%%%%% in manifold
from_R_to_01( A ) ->
    0.5 * (1.0 + math:atan(A) * 2 / math:pi() ).

% in F - target function, Man - manifold given as a function of singular cube
% out: {confirmed, [X1, X2, ...]} if local minimum found in given iterations, {presumably, [X1, X2, ...]} if not
minn_man(F, Man) ->
    MA = proplists:get_value( arity, erlang:fun_info( Man )),
    R1 = fun (X) -> from_R_to_01(X) end,

    OnRn_with_list = fun (L) ->
        OnMan = apply(Man, lists:map(R1, L)),
        apply(F, OnMan)
    end,
    OnRn = make_fun:nar(OnRn_with_list, MA),
    
    {Status, XY} = minn( OnRn ),
    {Status, apply(Man, lists:map(fun from_R_to_01/1, XY ) )}.


%%%%% in parallelotop
parallelotop( Pts ) ->
    Pn = length(Pts),
    Xn = length(hd(Pts)),
    ListTaking = fun( Xs ) -> 
        Pij = fun ( I, J ) -> lists:nth( J, lists:nth(I, Pts) ) end,
        Xi = fun ( I ) -> lists:nth( I, Xs ) end,
        [Pij(1, J) + lists:sum([Xi(I-1) * (Pij(I, J) - Pij(1, J)) || I <- lists:seq(2, Pn)]) || J <- lists:seq(1, Xn)]
    end,
    make_fun:nar(ListTaking, Pn-1).

% in: F - target function, Pts - parallelotop given by its basic points, 
%     eg [[0,0,0], [1,0,0], [0,1,0], [0,0,1]] would be 3-dimensional singular cube
% out: {confirmed, [X1, X2, ...]} if local minimum found in given iterations, {presumably, [X1, X2, ...]} if not
minn_para(F, Pts) ->
    minn_man( F, parallelotop( Pts )).


%%%%% in simplex
parallelotop_i (Sx, I) ->
    Xn = length(hd(Sx)),
    Switched = [lists:nth(I, Sx)] ++ lists:sublist(Sx, I-1) ++ lists:nthtail(I, Sx),
    X0 = fun (J) -> lists:nth(J, hd(Switched)) end,
    [[X0(J) + (lists:nth(J, X) - X0(J))/2.0 || J <- lists:seq(1, Xn) ] || X <- Switched].

f_pt_sx(F, L) ->
    lists:map(fun (X) -> minn_man(F, parallelotop(X)) end, L).

results_to_points(R) ->
    [P || {confirmed, P} <- R].

% in: F - target function, Sx - simplex given by its vertexes
% out: {confirmed, [X1, X2, ...]} if local minimum found in given iterations, {confirmed, []} if not
%      note, that minn_sx collects only confirmed minimums from its unity, therefore doesn't make presumptions
minn_sx(F, Sx) ->
    case length(Sx) of
        1 ->
            {confirmed, hd(Sx)};
        _ -> 
            Pts = [parallelotop_i(Sx, I) || I <- lists:seq(1, length(Sx))],
            Mins = f_pt_sx(F, Pts),
            {confirmed, minn_list(F, results_to_points( Mins ))}
    end.


%%%%% max - this section is just a bunch of declaration for corresponding minimum functions for your convenience

max1_list(F, L) ->
    min1_list( make_fun:inverse(F), L).

maxn_list(F, L) ->
    minn_list( make_fun:inverse(F), L).

max1(F) ->
    min1( make_fun:inverse(F)).

max1(F, I, E) ->
    min1( make_fun:inverse(F), I, E).

maxn(F) ->
    minn( make_fun:inverse(F)).

maxn(F, I, E) ->
    minn( make_fun:inverse(F), I, E).

maxn_man(F, M) ->
    minn_man( make_fun:inverse(F), M).

maxn_para(F, P) ->
    minn_para( make_fun:inverse(F), P).

maxn_sx(F, M) ->
    minn_sx( make_fun:inverse(F), M).


%%%%% nice display
nice_display(L) ->
    [erlang:display(X) || X <- L].


% test - simple test covering most of the functionality
t_f( X ) ->
    (X+5)*(X+5).

t_f2( X, Y ) ->
    (X-4)*(X-4) + (Y-3)*(Y-3).

man2(X, Y) ->
    [1.0 + 2*X - math:sin(Y), 0.5*math:sin(X) + 3*Y].

ptop() -> 
    [[4, 0], [1, 1], [0, 4]].

test() ->
    Sx = make_fun:from_text("[[1, 1], [0, 4], [4, 0]]", 0), % simplex
    Sx1 = make_fun:from_text("[[1, 1]]", 0),    % single point - 0-simplex
    Sxdg = make_fun:from_text("[[3, 1], [0, 4], [4, 0]]", 0),    % degenerative simplex
    Sxdge = make_fun:from_text("[[4, 0], [0, 4], [4, 0]]", 0),    % degenerative simplex with degenerative edge
    nice_display(["  Unconstrained:",
                  min1( fun t_f/1 ),
                  minn( fun t_f2/2 ),
                  "  On manifold:",
                  minn_man( fun t_f2/2, fun man2/2 ),
                  "  On paralelotop (including 'max..' test):",
                  maxn_para( make_fun:inverse( fun t_f2/2 ), ptop() ),
                  "  On simplexes:",
                  minn_sx( fun t_f2/2, Sx() ),
                  minn_sx( fun t_f2/2, Sx1() ),
                  minn_sx( fun t_f2/2, Sxdg() ),
                  minn_sx( fun t_f2/2, Sxdge() ),
                  "  the end."
     ]).

perf_test() ->
    {_, Go, _} = erlang:now(),
    N = length([ minn_sx( fun t_f2/2, ptop() ) || _ <- lists:seq(1, 100)]),
    {_, Now, _} = erlang:now(),
    io:format("~w in ~w~n", [N, Now-Go]).
